#' Select bulk samples and parent samples
#'
#' `select_sample_and_SNP()` select `highParent.GT`, `highParent.AD`, `highParent.GQ`, `lowParent.GT`, `lowParent.AD`, `lowParent.GQ`,
#' `highBulk.GT`, `highBulk.AD`, `highBulk.GQ`, `lowBulk.GT`, `lowBulk.AD` and `lowBulk.GQ` from a dataset generated from GATK' `VariantsToTable` module.
#' Effective SNP sites are selected according both parents and bulks information. Then bulk AD will be reassignment according parent genotype.
#'
#' @param data A data set generated by GATK's `VariantsToTable` module.
#' @param highP The name of high phenotype parent.
#' @param lowP The name of low phenotype parent.
#' @param highB The name of high phenotype bulk.
#' @param popType Population type, RIL or F2.
#' @param bulkSize A integer vector of two element, represent high bulk size and low bulk size.
#' @param lowB The name of low phenotype bulk.
#' @param minGQ Minimum GQ for parent.
#' @param refFreq SNP will be discarded if ERF Allelic Depth ratio is less than `refFreq` or more than `1-refFreq` in both bulks.
#'
#' @return A table SNP position, allelic depth.
#' @export
#'
#' @examples
#' library(easyQTLseq)
#' # Example with sample data from a GATK table.
#' file_path <- system.file("extdata", "subset.table.gz", package = "easyQTLseq")
#' # readr::read_tsv() has a faster speed than read.table() when reading a file.
#' data <- readr::read_tsv(file = file_path)
#' x <- select_sample_and_SNP(data = data, highP = "qY", lowP = "R3", highB = "Y", lowB = "R", popType = "F2", bulkSize = c(30, 30))
select_sample_and_SNP <- function(data, highP, lowP, highB, lowB, popType, bulkSize, minGQ = 10, refFreq = 0.3, chrLen){
  # 检查一些信息
  ## 检查refFreq
  if (refFreq >= 0.5 || refFreq <= 0) {
    stop("refFreq should > 0 and < 0.5. Exiting...\n")
  }
  ## 受否有chrLen，没有就生成
  if(missing(chrLen)) {
    chrLen <- data %>% group_by(CHROM) %>% summarise(Len = max(POS))
  }
  ## 当参考基因组是某个材料的
  #if (highP == "REF" || lowP == "REF") {
  #  data <- data %>% mutate(REF.GT = paste(REF, REF, sep = "/"))
  #}
  if ((!missing(highP) && highP=="REF") || (!missing(lowP) && lowP=="REF")) {
    data <- data %>% mutate(REF.GT = paste(REF, REF, sep = "/"))
  }
  # select bulk samples and parent samples (if exist)
  if (missing(highP) & missing(lowP)) {
    df1 <- data %>% select(
      CHROM, POS, REF, ALT,
      highBulk.GT = any_of(paste(highB, "GT", sep = ".")),
      lowBulk.GT = any_of(paste(lowB, "GT", sep = ".")),
      highBulk.AD = any_of(paste(highB, "AD", sep = ".")),
      lowBulk.AD = any_of(paste(lowB, "AD", sep = ".")),
      highBulk.GQ = any_of(paste(highB, "GQ", sep = ".")),
      lowBulk.GQ = any_of(paste(lowB, "GQ", sep = "."))
    ) %>% na.omit()
  } else if (missing(highP)) {
    df1 <- data %>% select(
      CHROM, POS, REF, ALT,
      lowParent.GT = any_of(paste(lowP, "GT", sep = ".")),
      highBulk.GT = any_of(paste(highB, "GT", sep = ".")),
      lowBulk.GT = any_of(paste(lowB, "GT", sep = ".")),
      lowParent.AD = any_of(paste(lowP, "AD", sep = ".")),
      highBulk.AD = any_of(paste(highB, "AD", sep = ".")),
      lowBulk.AD = any_of(paste(lowB, "AD", sep = ".")),
      lowParent.GQ = any_of(paste(lowP, "GQ", sep = ".")),
      highBulk.GQ = any_of(paste(highB, "GQ", sep = ".")),
      lowBulk.GQ = any_of(paste(lowB, "GQ", sep = "."))
    ) %>% na.omit()
  } else if (missing(lowP)) {
    df1 <- data %>% select(
      CHROM, POS, REF, ALT,
      highParent.GT = any_of(paste(highP, "GT", sep = ".")),
      highBulk.GT = any_of(paste(highB, "GT", sep = ".")),
      lowBulk.GT = any_of(paste(lowB, "GT", sep = ".")),
      highParent.AD = any_of(paste(highP, "AD", sep = ".")),
      highBulk.AD = any_of(paste(highB, "AD", sep = ".")),
      lowBulk.AD = any_of(paste(lowB, "AD", sep = ".")),
      highParent.GQ = any_of(paste(highP, "GQ", sep = ".")),
      highBulk.GQ = any_of(paste(highB, "GQ", sep = ".")),
      lowBulk.GQ = any_of(paste(lowB, "GQ", sep = "."))
    ) %>% na.omit()
  } else {
    df1 <- data %>% select(
      CHROM, POS, REF, ALT,
      highParent.GT = any_of(paste(highP, "GT", sep = ".")),
      lowParent.GT = any_of(paste(lowP, "GT", sep = ".")),
      highBulk.GT = any_of(paste(highB, "GT", sep = ".")),
      lowBulk.GT = any_of(paste(lowB, "GT", sep = ".")),
      highParent.AD = any_of(paste(highP, "AD", sep = ".")),
      lowParent.AD = any_of(paste(lowP, "AD", sep = ".")),
      highBulk.AD = any_of(paste(highB, "AD", sep = ".")),
      lowBulk.AD = any_of(paste(lowB, "AD", sep = ".")),
      highParent.GQ = any_of(paste(highP, "GQ", sep = ".")),
      lowParent.GQ = any_of(paste(lowP, "GQ", sep = ".")),
      highBulk.GQ = any_of(paste(highB, "GQ", sep = ".")),
      lowBulk.GQ = any_of(paste(lowB, "GQ", sep = "."))
    ) %>% na.omit()
  }
  # 如果有些GT是“A|A”形式的，将“|”换成“/”
  df1 <- df1 %>%
    mutate(across(ends_with(".GT"), ~ str_replace(.x, "\\|", "/")))
  # filter according parent GT and GQ
  ## filter GQ
  df2 <- df1
  if ("highParent.GQ" %in% colnames(df2)) {
    df2 <- df2 %>% filter(highParent.GQ > minGQ)
  }
  if ("lowParent.GQ" %in% colnames(df2)) {
    df2 <- df2 %>% filter(lowParent.GQ > minGQ)
  }
  ## fitler GT
  if ("highParent.GT" %in% colnames(df2)) {
    df2 <- df2 %>% filter(highParent.GT == paste(REF, REF, sep = "/") | highParent.GT == paste(ALT, ALT, sep = "/"))
  }
  if ("lowParent.GT" %in% colnames(df2)) {
    df2 <- df2 %>% filter(lowParent.GT == paste(REF, REF, sep = "/") | lowParent.GT == paste(ALT, ALT, sep = "/"))
  }
  if (all(c("highParent.GT", "lowParent.GT") %in% colnames(df2))) {
    df2 <- df2 %>% filter(highParent.GT != lowParent.GT)
  }

  # 将highBulk和lowBulk的AD分开，并计算DP值
  df3 <- df2 %>%
    separate(highBulk.AD, c("highBulk.AD_0", "highBulk.AD_1"), sep = ",", convert = TRUE) %>%
    separate(lowBulk.AD, c("lowBulk.AD_0", "lowBulk.AD_1"), sep = ",", convert = TRUE) %>%
    mutate(
      HB.DP = highBulk.AD_0 + highBulk.AD_1,
      LB.DP = lowBulk.AD_0 + lowBulk.AD_1
    )
  if ("highParent.AD" %in% colnames(df3)) {
    df3 <- df3 %>%
      separate(highParent.AD, c("highParent.AD_0", "highParent.AD_1"), sep = ",", convert = TRUE) %>%
      mutate(
        HP.DP = highParent.AD_0 + highParent.AD_1
      )
  }
  if ("lowParent.AD" %in% colnames(df3)) {
    df3 <- df3 %>%
      separate(lowParent.AD, c("lowParent.AD_0", "lowParent.AD_1"), sep = ",", convert = TRUE) %>%
      mutate(
        LP.DP = lowParent.AD_0 + lowParent.AD_1
      )
  }

  # 过滤HB index和LB index同时高于某值或同时低于某值的位点
  df4 <- df3 %>%
    filter(!(((highBulk.AD_0 / HB.DP < refFreq) &
                (lowBulk.AD_0 / LB.DP < refFreq)) |
               ((highBulk.AD_0 / HB.DP > 1-refFreq) &
                  (lowBulk.AD_0 / LB.DP > 1-refFreq)))
    )

  # 根据亲本、混池基因型，将混池的AD_0和AD_1重新分配
  if ("highParent.GT" %in% colnames(df4)) {
    df5 <- df4 %>%
      mutate(
        HB.HP.AD = if_else(highParent.GT == paste(REF, REF, sep = "/"),
                           highBulk.AD_0, highBulk.AD_1),
        HB.LP.AD = if_else(highParent.GT == paste(REF, REF, sep = "/"),
                           highBulk.AD_1, highBulk.AD_0),
        LB.HP.AD = if_else(highParent.GT == paste(REF, REF, sep = "/"),
                           lowBulk.AD_0, lowBulk.AD_1),
        LB.LP.AD = if_else(highParent.GT == paste(REF, REF, sep = "/"),
                           lowBulk.AD_1, lowBulk.AD_0)
      ) %>%
      select(
        CHROM, POS, REF, ALT,
        any_of(c("HP.DP", "LP.DP")),
        HB.HP.AD, HB.LP.AD, HB.DP,
        LB.HP.AD, LB.LP.AD, LB.DP
      )
  } else if ("lowParent.GT" %in% colnames(df4)) {
    df5 <- df4 %>%
      mutate(
        HB.HP.AD = if_else(lowParent.GT == paste(ALT, ALT, sep = "/"),
                           highBulk.AD_0, highBulk.AD_1),
        HB.LP.AD = if_else(lowParent.GT == paste(ALT, ALT, sep = "/"),
                           highBulk.AD_1, highBulk.AD_0),
        LB.HP.AD = if_else(lowParent.GT == paste(ALT, ALT, sep = "/"),
                           lowBulk.AD_0, lowBulk.AD_1),
        LB.LP.AD = if_else(lowParent.GT == paste(ALT, ALT, sep = "/"),
                           lowBulk.AD_1, lowBulk.AD_0)
      ) %>%
      select(
        CHROM, POS, REF, ALT,
        any_of(c("HP.DP", "LP.DP")),
        HB.HP.AD, HB.LP.AD, HB.DP,
        LB.HP.AD, LB.LP.AD, LB.DP
      )
  } else {
    df5 <- df4 %>%
      select(CHROM, POS, REF, ALT,
             HB.REF.AD = highBulk.AD_0, HB.ALT.AD = highBulk.AD_1, HB.DP,
             LB.REF.AD = lowBulk.AD_0, LB.ALT.AD = lowBulk.AD_1, LB.DP)
  }

  if (!missing(highP) && !missing(lowP)) {
    QTLseq <- structure(
      list(data = df5, highP = highP, lowP = lowP, highB = highB, lowB = lowB,
           popType = popType, bulkSize = bulkSize, slidwin = data.frame(), chrLen = chrLen),
      class = c("QTLseq", "WithParent", "BothParent")
    )
  } else if (!missing(highP)) {
    QTLseq <- structure(
      list(data = df5, highP = highP, lowP = "lowParent", highB = highB, lowB = lowB,
           popType = popType, bulkSize = bulkSize, slidwin = data.frame(), chrLen = chrLen),
      class = c("QTLseq", "WithParent", "HighParent")
    )
  } else if (!missing(lowP)) {
    QTLseq <- structure(
      list(data = df5, highP = "highParent", lowP = lowP, highB = highB, lowB = lowB,
           popType = popType, bulkSize = bulkSize, slidwin = data.frame(), chrLen = chrLen),
      class = c("QTLseq", "WithParent", "LowParent")
    )
  } else {
    QTLseq <- structure(
      list(data = df5, highB = highB, lowB = lowB,
           popType = popType, bulkSize = bulkSize, slidwin = data.frame(), chrLen = chrLen),
      class = c("QTLseq", "WithoutParent")
    )
  }
  return(QTLseq)
}
